#!/bin/csh -f
set echo_style=both

# 1. first set of default values is below
# 2. second set is read from INPUT/DFIT/defaults
# 3. third set are commandline arguments
# RULE: the later set values OVERRIDE: cmdline overrides defaults from file, and defaults from file override defaults set in this script

# dfit: diffraction/relaxation
set commandline="$argv"

set SRCDIR="$DECAHOME/SRC"
set EXEC="$DECAHOME/bin/dfit_"
set PHELP="awk -f $SRCDIR/pophelp.awk $SRCDIR/dfit.help"
set LSOPT="$DECAHOME/scripts/lsopt dfit"
set CROMERNAME="$DECAHOME/DiffData/cromer.scf"

#### DEFAULT TILING
set FNAME="tiling.inp"

#### DIFFRACTION DATA
set dpathdef="$DECAHOME/DiffData"
set dpath=$dpathdef
#set dname="alcocu"
set dname="alnico"
set dname="alnico_ex002"
#set dname="anc.c2s"
#set dname="anc.c2s_ex002"
#set dname="anc-app.e2s.hkl"
#set dname="anc-app-lift.e2s.hkl"
#set dname="anc-app2deca.hkl"

#### INITIALIZATION PARAMETERS
# global

set d04inp=0 # this option should be set by d04 only

set ORB_LIST=()  # N O_1 ... O_N: explicit list of orbits to refine/relax
# relaxation
set RCPLUS=1.0  # [Angs] add to rcut for construction of neighbor tables
# MONTE CARLO
# note: for MC, I used lambda=0.2
set lambda=0.0001 # good for ueff=0 and Steurers old AlNiCo data
set etarg="-0.35"
set equil=0.5
set sigdr=0.1 # in Angstrom
set cueff=0   # bodn-weighting for energy
set DMCS=0 # MC
set CMCS=0 # MC
set hloops=1
set loops=1
set dmcupmode=2 # dw update mode for dmc
set efitmode=2
# diffraction
set INFO=0
set ITER=1       #  dfit iterations within a single relaxation move
set ALAMU=10.0
set DCHI=0.1

set et1="-0.4"
set et2="-0.3"
set la1="0.0001"
set la2="0.0001"
set be1=10
set be2=50
set fita=()
set nfit=0
set apsig=(0 0)

set udw=()
set nsigdw=0
set sigdw=()  # for DW-ave datapoints
# LIMITS
set dwlim=(-10 100)
set adwlim=(-10 100)
set dwplim=(-10 100)
set aclim=(0 100.0)
set dwglo=0

set DW_RELOAD=0   #  1->reload DW factors from INPUT/DFIT/ini
set ANISO=0       #  refine anisotropic DW's?
set DWPERP=0      #  [0=none 1=global]
set ACONST=0      #  refine additive constant F_obs=sqrt(I_calc+const)
set QPLIM=5.0     #  Q-perp limit on the expe data
set APLIM=100     #  elimination of bad approximant peaks
set aphkl="none"
set SIGMALIM=5.0  #  limit on F_obs/sigma
set FLIM=0.0  #  limit on F_obs/sigma
set ONLYEVEN=0   # use only fundamental datapoints (no superstructure)

######## set defaults from INPUT/DFIT/defaults ########
set def="INPUT/DFIT/defaults"
set argv_def=""
if ( -e $def ) then
    set argv_def=`awk 'NF>0{sub(/^ */,"");a=index($0,"#");c=$0;if(a>1){split($0,b,"#");c=b[1]}if(a!=1){cmd=cmd" -"c}}END{printf cmd}' $def`
#    echo $argv_def
endif
########################### OPTIONS ###############################
set TEST=0
set HELP=0
set argv_save=($argv)
set argv=($argv_def)
set REREAD=1
startoptions:
#echo "reread $REREAD"
#echo "$argv"
while( "$*" != "" )
#    echo $1
    if ( -e $1 ) then # config name is not option
	set FNAME=$1
	shift
    endif
    if ( $# == 0 ) then
	goto skip_options
    endif
    switch ( $1 )
	case -data:
	    shift
	    set dname=$1
	breaksw
	case -path:
	    shift
	    set dpath=$1
	breaksw
	case -test:
	    set TEST=1
	breaksw
	case -d04inp:
	    set d04inp=1
	    set TEST=1
	breaksw
	case -dwglo:
	    @ dwglo = 1 - $dwglo
	breaksw
	case -rdw:
	    set DW_RELOAD = 1
	    set pdump=0
	breaksw
	case -para:
	    shift
	    set DW_RELOAD=2
	    set pdumpname=$1
	    if (! -e $pdumpname ) then
		echo "-pdump: file $pdumpname does not exist"
		exit
	    endif
	    echo "Reloading DW factors from $pdumpname"
	breaksw
	case -dwperp:
	    set DWPERP = 1
	breaksw
	case -dwperp2:
	    set DWPERP = 2
	breaksw
	case -dwlim:
	    shift
	    set dwlim=($1) # lower bound
	    shift
	    set dwlim=($dwlim $1) # upper bound
	breaksw
	case -dmcupmode:
	    shift
	    set dmcupmode=$1
	breaksw
	case -adwlim:
	    shift
	    set adwlim=($1) # lower bound
	    shift
	    set adwlim=($adwlim $1) # upper bound
	breaksw
	case -dwplim:
	    shift
	    set dwplim=($1) # lower bound
	    shift
	    set dwplim=($dwplim $1) # upper bound
	breaksw
	case -aclim:
	    shift
	    set aclim=($1) # lower bound
	    shift
	    set aclim=($aclim $1) # upper bound
	breaksw
	case -udw:
	    shift
	    set udw=($1)
	breaksw
	case -udwch:
	    shift
	    set sigdw=($sigdw $1)  # chemistry
	    shift
	    set sigdw=($sigdw $1)  # value
	    @ nsigdw++
	breaksw
	case -dmc:
	    shift
	    set DMCS=$1
	breaksw
	case -cmc:
	    shift
	    set CMCS=$1
	breaksw
	case -hl:
	    shift
	    set hloops=$1
	breaksw
	case -lambda:
	    shift
	    set la1=$1
	    set la2=$1
	breaksw
	case -beta:
	    shift
	    set be1=$1
	    set be2=$1
	breaksw
	case -sigdr:
	    shift
	    set sigdr=$1
	breaksw
	case -ueff:
	    shift
	    set cueff=$1
	breaksw
	case -fit:
	    shift
	    set f1=$1
	    shift
	    set f2=$1
	    set fita=($fita $f1 $f2)
	    @ nfit++
	breaksw
	case -nofit:
	    set fita=()
	    set nfit=0
	breaksw
	case -lambdas:
	    shift
	    set la1=$1
	    shift
	    set la2=$1
	breaksw
	case -betas:
	    shift
	    set be1=$1
	    shift
	    set be2=$1
	breaksw
	case -etargs:
	    shift
	    set et1=$1
	    shift
	    set et2=$1
	breaksw
	case -etarg:
	    shift
	    set et1=$1
	    set et2=$1
	breaksw
	case -loops:
	    shift
	    set loops=$1
	breaksw
	case -rcplus:
	    shift
	    set RCPLUS=$1
	breaksw
	case -iter:
	    shift
	    set ITER=$1
	breaksw
        case -alamu:
            shift
	    set ALAMU=$1
        breaksw
	case -dchi:
	    shift
	    set DCHI=$1
	breaksw
	case -sglim:
	    shift
	    set SIGMALIM=$1
	breaksw
	case -flim:
	    shift
	    set FLIM=$1
	breaksw
	case -qplim:
	    shift
	    set QPLIM=$1
	breaksw
	case -apsig:
	    shift
	    set apsig=(1 $1)
	breaksw
	case -aplim:
	    shift
	    set aplim=$1
 	    shift
	    set aphkl=$1
	breaksw
	case -aconst:
	    set ACONST = 1
	breaksw
	case -aniso:
	    shift
	    set ANISO = $1
	    if( ($ANISO < 0) || ($ANISO > 3))then
		echo "option -aniso: ANISO must be intreger >=0 and <=3"
		exit
	    endif
	breaksw
	case -info:
	    set INFO = 1
	breaksw
	case -olist:
	    shift
	    if ( -e $1 ) then 
		set ORB_LIST=($ORB_LIST `head -1 $1`)
	    else
		echo "Option -olist [FILE]: file $1 with list of esg values not found."
		exit
	    endif
	breaksw
	case -h:
	    goto helpmessage
	breaksw
        case -help:
	    helpmessage:
            shift
            if ( $# == 0 ) then
		$LSOPT
		exit
            endif
	    $PHELP $1
            exit
        breaksw
	default:
	    echo "Unrecognized option $1"
	    $LSOPT
	    exit
	breaksw
    endsw
    shift
end
if ( $REREAD == 1 ) then
    set argv=`echo $argv_save`
    set REREAD=0
    if ( $# > 0 ) goto startoptions
endif

skip_options:
#### END OPTIONS
if !( -e $FNAME ) then 
  echo "dfit: tiling/config input file $FNAME does not exist." > /dev/stderr
  exit
endif
# check if the input file is type copa.dfit
if ( $DW_RELOAD != 1 ) then
    set parainp=$FNAME # default, for $DW_RELOAD==0
    if( $DW_RELOAD == 2 ) then
	set parainp=$pdumpname
    endif
    set has_dfit_pars=`awk '$1=="start_dfit_params"{a=1}$1=="end_dfit_params"{b=1}END{if(a&&b){printf "1"}else{printf "0"}}' $parainp`
#    echo "$has_dfit_pars $parainp"
    set ppardump="INPUT/DFIT/ppar.dump"
    if ( $has_dfit_pars == 1 ) then
	if($DW_RELOAD == 0) then 
	    echo "$FNAME has dfit section: reloading DW paramters (disable by -rdw -1)"
	    $DECAHOME/scripts/printdfpars $parainp > $ppardump
	    set DW_RELOAD=3
	endif
    else if ( $DW_RELOAD == 2 ) then
	echo "Warning: could not find dfit parameters section in $parainp: skipping reload"
	set DW_RELOAD=0
    endif
endif
# diffraction stuff
if !( -e $CROMERNAME ) then
  echo "Missing file ENV(DECAHOME)/DiffData/cromer.scf." > /dev/stderr
  exit
endif
if (( $DW_RELOAD == 1 ) && !( -e INPUT/DFIT/ini )) then
  echo "LOAD_DW=1: missing file INPUT/DFIT/ini" > /dev/stderr
  exit
endif
if (( $DW_RELOAD == 1 ) && ( $INFO == 1 )) then
    echo "Reloading DW parameters from INPUT/DFIT/ppar.ini" > /dev/stderr
endif
set DIFFDATA=$dpath"/"$dname
if !( -r $DIFFDATA ) then 
    echo "Diffraction-DataFile $DIFFDATA does not exist." > /dev/stderr
    echo "Full path to the default diff-data dir is $dpathdef." > /dev/stderr
    exit
endif
set nelim=0
# The approximant data contain word APPROXIMANT on the first line ->
set APPROX="`grep -c APPROX $DIFFDATA`"   #returns 1 for approximant data
# Check diffraction data file
if( $APPROX == "0" ) then
    #TEST APP DIFFRACTION DATA FILE FIRST, eliminate multiply defined datapts
    $DECAHOME/scripts/data2app -mode 2 -data $dname $FNAME -test
    awk 'FNR>4{p[$5]++;p[$6]++}END{for(i in p){print i}}' test_hkl.dat | sort -n > INPUT/DFIT/hkl.elim
    if( -e $aphkl) then
	awk 'BEGIN{x=ARGV[2];delete ARGV[2]; ARGC--}$4/$2>x{print $5}' $aphkl $aplim >> INPUT/DFIT/hkl.elim
    endif
    set nout=`wc -l INPUT/DFIT/hkl.elim`
    if( $nout[1]>0 ) then
	echo "$nout[1] bad approximant datapoints will be eliminated"
    endif
    mv test_hkl.dat INPUT/DFIT
    $DECAHOME/scripts/data2app -mode 1 -data $dname $FNAME > INPUT/DFIT/dat.hkl  #generates dat.hkl
    if ($INFO == 1) then 
	echo "DiffractionData: $dname"  > /dev/stderr
    endif
    set nelim=$nout[1]
else
  if( $ONLYEVEN == 1 ) then 
    echo "Using EVEN approximant data from the dataset $dname"  > /dev/stderr
    awk 'FNR<3{print}FNR>2&&$5%2==0{print}' $DIFFDATA > INPUT/DFIT/dat.hkl
  else
    echo "Copying approximant data $dname onto INPUT/DFIT/dat.hkl" > /dev/stderr
    cp $DIFFDATA INPUT/DFIT/dat.hkl
  endif
endif

#AUTOMATIC SETING OF INPUT_TYPE AND CONSISTENCY CHECKS:
if !( -e $SRCDIR/tico.awk ) then
    echo "dfit: File $SRCDIR/tico.awk is required." > /dev/stderr
    exit
endif
set FNOK=`awk -f $SRCDIR/tico.awk -s $FNAME`
set TINP=$FNOK[1]
if ( $TINP > 1 ) then
    echo "TINP>1: $FNAME is neither TILING nor CONFIG type." > /dev/stderr
    exit
endif
set ERR=$FNOK[2]
set ERRPO=3
if ( $TINP > 0 ) then
    set ERRPO=4
endif
set ERRMSG=`echo $FNOK[$ERRPO-]`
set TMSG=`echo "TILING CONFIG FLIP_REP"`
if ( $ERR > 0 ) then
    @ ind = $TINP + 1
    echo $FNAME" ("$TMSG[$ind]" type): "$ERRMSG > /dev/stderr
    exit
endif
set out="INPUT/dfit.inp"
if ($TEST == 1)then 
    set out="/dev/stdout"
    echo "-------------------- start dfit.inp "
endif
if($loops > 1 ) then
    set LAMBDAS=(`awk -f $DECAHOME/SRC/dfit_setesg.awk $loops $la1 $la2`)
    set ETARGS=(`awk -f $DECAHOME/SRC/dfit_settem.awk $loops $et1 $et2`)
    set BETAS=(`awk -f $DECAHOME/SRC/dfit_settem.awk $loops $be1 $be2`)
else
    set LAMBDAS=($la1)
    set ETARGS=($et1)
    set BETAS=($be1)
endif
#launch!
if ( $nfit == 0 && $DMCS == 0 && $CMCS == 0 ) set loops=0
echo "$nfit  $fita # fit # fit actions " > $out
echo "$loops $hloops # loops # hloops #" >> $out
echo "$LAMBDAS # lambdas #" >> $out
echo "$ETARGS  # etargs #" >> $out
echo "$BETAS  # betas #" >> $out
echo "$DMCS $CMCS $dmcupmode # dmcs cmcs dmcdwupmode" >> $out
echo "$sigdr  # sigdr # " >> $out
echo "$cueff # cueff" >> $out
echo "$dwglo  # dwglo # fit elemental dw factors instead orbit DW factors" >> $out
echo "$ANISO  # aniso # int ANISO: 1 = refine anisotropic DW factors" >> $out
echo "$DWPERP  # dwperp # int DWPERP: 1 = refine one global DW-perp factor" >> $out
echo "$ACONST  # aconst # int ACONST: 1 = refine additive constant F_obs=sqrt(I_calc+const)" >> $out
echo "$SIGMALIM  # sglim # real SIGMALIM: limit on F_obs/sigma" >> $out
echo "$FLIM  # flim # real FLIM: limit on minimal F_obs" >> $out
echo "$QPLIM  # qplim # real QPLIM: Q-perp limit on the expe data" >> $out
echo "$dwlim  # upper/lower limits for deca/global DW factors" >> $out
echo "$adwlim  # upper/lower limits for aniso/Z DW factors" >> $out
echo "$dwplim  # upper/lower limits for perp DW factor" >> $out
echo "$aclim  # upper/lower limits background const subtraction" >> $out
echo "$#udw $udw  # set accuracy with which DWs are uniform for all DWs">>$out
echo "$nsigdw $sigdw  # CHEM accuracy for a specie" >> $out
echo "$ITER  $ALAMU $DCHI # iter alamu dchi # int ITER: real ALAMU, DCHI : MRQ setup " >> $out
echo "$#ORB_LIST $ORB_LIST # olist # ORB_LIST: list of refined orbits; none->ALL" >> $out
echo "$DW_RELOAD  # rdw # int DW_RELOAD: 1 = reload DW factors from INPUT/DFIT/ini" >> $out
echo "$nelim  # eliminate ambiguous datapts (autoset) from INPUT/DFIT/hkl.elim" >> $out
echo "$apsig  # add app-sigma to the error bars" >> $out
echo "$RCPLUS  # rcplus # real RCLPLUS [A]: increase cutoff for construction of neighbor tables" >> $out
echo "$CROMERNAME" >> $out
echo "$INFO  # info # int INFO: set verbose" >> $out
#LAUNCH (when not test)
if (! -e DFIT) mkdir DFIT
if($TEST == 1)then
    if ( $d04inp == 0 ) then 
	echo "-------------------- pipe to exec: "
	echo $TINP\\n$FNAME
    endif
else # LAUNCH!
    echo "dfit "$commandline >> INPUT/dfit.cmdline
    echo $TINP\\n$FNAME | $EXEC
endif
##

